"""
三维克里金插值
    1.data1的输入路径需要为 poin_and_depth_data.xlsx的路径
    2.input_points的输入路径需要为GridB_change.xlsx的路径
    
"""



import pandas as pd
import numpy as np
from skgstat import Variogram, OrdinaryKriging
from scipy.spatial import cKDTree
from concurrent.futures import ThreadPoolExecutor, as_completed
from tqdm import tqdm
# 读取第一份文件数据
data1 = pd.read_excel(r"D:\360MoveData\Users\10758\Desktop\yuanzi\yuanzi_code\griddatatrans-master-83f28c8ab3938ed6aa1007f0e027b365e03c0950\yuanzi\q1_result\poin_and_depth_data.xlsx")
coordinates1 = data1[['x', 'y', 'depth']].values.astype(float)  # 转换为float类型
young_modulus = data1['young_modulus'].values.astype(float)
poisson_ratio = data1['poisson_ratio'].values.astype(float)
# 创建KD树
tree = cKDTree(coordinates1)
# 读取第二份文件数据
input_points = pd.read_excel(r"D:\360MoveData\Users\10758\Desktop\yuanzi\yuanzi_code\griddatatrans-master-83f28c8ab3938ed6aa1007f0e027b365e03c0950\yuanzi\q1_result\GridB_change.xlsx")
input_coordinates = input_points[['x', 'y', 'depth']].values.astype(float)
# 定义处理单个点的函数
def process_point(input_point, tree, coordinates1, values):
    # 使用 KD 树查找与输入点距离在 r 以内的所有点
    indices = tree.query_ball_point(input_point, r=80)  # 查询半径为 r 的点的索引
    if not indices:  # 如果没有找到任何邻近点，则返回 NaN
        return np.nan
    # 限制邻居点数量z
    min_neighbors = 50  # 最小邻近点数
    max_neighbors = 100  # 最大邻近点数
    if len(indices) > max_neighbors:
        distances, indices = tree.query(input_point, k=max_neighbors)   
    # 根据找到的邻近点提取杨氏模量或泊松比
    neighbor_values = values[indices]
    neighbor_coordinates = coordinates1[indices] 
    # 检查是否有足够的邻居点
    if len(neighbor_values) < min_neighbors:  # 确保至少有 30 个邻居点
        print(f"Warning: for point {input_point}, not enough neighbors were found within the range.")
        return np.nan  # 或者返回一个默认值
    # 创建变差函数模型
    try:
        V = Variogram(
            neighbor_coordinates,
            neighbor_values, 
            model='spherical',  # 可以尝试其他模型，如 "spherical",'exponential', 'gaussian' 等
            normalize=False,  # 不进行归一化处理
            n_lags=1,  # 增加滞后数 #(细化变差函数的分辨率)
            maxlag=None, # None
            verbose=True  # 打印详细信息，帮助调试
        ) 
        
        # 增加最大迭代次数
        V.fit(force=True, bounds=None, p0=None, maxfev=50000000)  # 增加 maxfev       
        # 构建 Kriging 模型
        ok = OrdinaryKriging(V,min_points=1)      
        # 进行插值预测
        try:
            predicted_value = ok.transform(input_point.reshape(1, -1))
        except Exception as e:
            # 记录错误并设置预测值为 NaN
            print(f"Error at point {input_point}: {e}")
            predicted_value = np.nan
    except RuntimeError as re:
        print(f"RuntimeError creating variogram for point {input_point}: {re}")
        return np.nan  # 变差函数创建失败，返回 NaN
    except Exception as e:
        print(f"Error creating variogram for point {input_point}: {e}")
        return np.nan  # 其他异常，返回 NaN
    return predicted_value[0] if isinstance(predicted_value, np.ndarray) else predicted_value
# 使用 ThreadPoolExecutor 来并行处理每个点
max_workers = 10 # 你可以根据你的硬件条件调整这个数值
with ThreadPoolExecutor(max_workers=max_workers) as executor:
    # 提交所有任务
    futures_young = [executor.submit(process_point, point, tree, coordinates1, young_modulus) for point in input_coordinates]
    futures_poisson = [executor.submit(process_point, point, tree, coordinates1, poisson_ratio) for point in input_coordinates]  
    # 收集结果
    predicted_youngs = []
    for future in tqdm(as_completed(futures_young), total=len(futures_young), desc="Processing Young's modulus points"):
        predicted_youngs.append(future.result())
    predicted_poissons = []
    for future in tqdm(as_completed(futures_poisson), total=len(futures_poisson), desc="Processing Poisson's ratio points"):
        predicted_poissons.append(future.result())
# 将结果保存到 DataFrame 并导出
results = pd.DataFrame({
    'X': input_coordinates[:, 0],
    'Y': input_coordinates[:, 1],
    'Depth': input_coordinates[:, 2],
    'Predicted Young Modulus': predicted_youngs,
    'Predicted Poisson Ratio': predicted_poissons
})
results.to_excel(r"D:\360MoveData\Users\10758\Desktop\yuanzi\yuanzi_code\griddatatrans-master-83f28c8ab3938ed6aa1007f0e027b365e03c0950\q1_code\allresult_spherical_minpoint1_maxlagNone_nlags1.xlsx", index=False)
print("所有预测完成，结果已保存。")